Mouse lung BLM + AB RNA-seq analysis

Intro

This is the RNA-seq differential expression analysis for a prevention experiment using the X203 antibody during 3 weeks of bleomycin treatment in mice. Lungs are harvested, the RNA is extracted and sequenced. The 3 conditions are Saline (no bleomycin), IgG (bleomycin + IgG, control) and X203 (bleomycin + X203, treatment).

Library loading

libs = c("DESeq2", "sva", "edgeR", "tidyverse", "biomaRt", "beeswarm", "gplots", "pheatmap", "viridis", "treemapify", "ggplot2", "scales", "colorspace", "MSigDB", "fgsea")
suppressMessages({lapply(libs, require, character.only = TRUE)})
source("helper_functions.R")
while (dev.cur()>1) dev.off()

Data loading

The gene-level count table and the sample names table are loaded. Counts are formatted and the coldata object is created.

#Load the counttable and sample names files
counttable <- read.delim("Lung_mouse_all.counts", header = T, comment.char = '#', stringsAsFactors = F)
samplesheet <- read.delim("samplesheet.txt", header = T, stringsAsFactors = F)
colnames(counttable)[7:18] <- gsub(colnames(counttable)[7:18], pattern = ".Aligned.sortedByCoord.out.bam", replacement = "")
counttable = counttable[,c(1,6,7:18)]

#Row names as sample ID in samplesheet
rownames(samplesheet) <- samplesheet$Sample_ID

#Change colnames of counttable to Description from samplesheet
colnames(counttable)[3:14] = samplesheet[colnames(counttable)[3:14],7]

#Use geneid as rownames in counttable table
rownames(counttable) = counttable$Geneid
counttable$Geneid = NULL

#Make coldata
coldata = as.data.frame(colnames(counttable)[2:13])
coldata$treatment = rep(c("Saline", "IgG", "X203"), 4)
coldata$batch = rep(c("1","2","3","4"), each = 3)
colnames(coldata)[1] = "sample"
rownames(coldata) = coldata$sample
coldata$treatmentColor = rep(c("gray", "purple", "salmon"), 4)
colnames(counttable)[2:13] = paste(coldata$treatment, coldata$batch, sep = "_")

#Obtain gene symbols from biomart
ensembl = useMart("ensembl", dataset = "mmusculus_gene_ensembl")
mapping <- getBM(attributes = c("ensembl_gene_id", "mgi_symbol"), mart = ensembl)
mapping <- mapping[-which(duplicated(mapping[,1])),]
rownames(mapping) <- mapping[,1]
mapping2 <- mapping[rownames(counttable),]

#Get the orthologs between human and mouse
mart <- useDataset("mmusculus_gene_ensembl", mart=useMart("ensembl"))
bm <- getBM(attributes=c("ensembl_gene_id", "hsapiens_homolog_associated_gene_name"), mart=mart) %>%
  distinct() %>%
  as_tibble() %>%
  na_if("") %>% 
  na.omit()
bm = as.data.frame(bm)

Quality control

Number of reads

par(mar = c(5,5,2,2))
barplot(colSums(counttable[,2:13]), las = 2, border = NA, col = "slateblue", main = "Gene-level count sample size", ylab = "counts")

We have plenty of counttable to perform direffential expression accurately.

PCA

genetable <- data.frame(gene.id=rownames(counttable))
z <- DGEList(counts=counttable[,2:13], samples=coldata, genes=genetable)
my_mod = model.matrix(~ z$samples$treatment)
All_norm <- calcNormFactors(z, method="upperquartile")
my_data = cpm(All_norm, log=TRUE, prior.count=2)
plot_PCA_2(my_data, 1, 2, coldata, coldata$treatmentColor, symbol = 16, main = "PCA of log2(CPM)", lab1 = coldata$treatment, lab2 = coldata$batch)

The PCA shows that the first principal component separates properly the 3 treatments, and that the X203 treatment is intermediate between Saline and IgG. However, the position of X203_3 suggests that it is an outlier. Since we have at least 3 samples per treatment, we can remove it. However we first check Cook’s distance distributions:

rownames(coldata) <- colnames(counttable[,2:13])
ddsmat_all = DESeqDataSetFromMatrix(countData = counttable[,2:13], colData = coldata, design = ~ treatment)

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

#Estimate size factors: necessary to generate a normalized count table
ddsmat_all = estimateSizeFactors(ddsmat_all)

#Remove genes with 0 counttable in all samples to speed computation up
ddsmat_all_ne = ddsmat_all[rowSums(counts(ddsmat_all)) > 1, ]

#Differential expression
dds = DESeq(ddsmat_all_ne)

## using pre-existing size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

boxplot(log10(assays(dds)[["cooks"]]), border = c(rep("black", 8), "red", rep("black", 3)))

Although it does not seem to be extremely out of line, we remove it.

Outlier removal

counttable_2 = counttable[,-10]

#Make coldata
coldata_2 = as.data.frame(colnames(counttable_2)[2:12])
coldata_2$treatment = c("Saline", "IgG", "X203", "Saline", "IgG", "X203", "Saline", "IgG", "Saline", "IgG", "X203")
colnames(coldata_2)[1] = "sample"
rownames(coldata_2) = coldata_2$sample
coldata_2$treatmentColor = c("gray", "purple", "salmon", "gray", "purple", "salmon", "gray", "purple", "gray", "purple", "salmon")



genetable <- data.frame(gene.id=rownames(counttable_2))
z <- DGEList(counts=counttable_2[,2:12], samples=coldata_2, genes=genetable)
my_mod = model.matrix(~ z$samples$treatment)
All_norm <- calcNormFactors(z, method="upperquartile")
my_data = cpm(All_norm, log=TRUE, prior.count=2)

 plot_PCA_2(my_data, 1, 2, coldata_2, coldata_2$treatmentColor, symbol = 16, main = "PCA of log2(CPM)", lab1 = coldata_2$treatment, lab2 = coldata_2$batch)

Without the outlier, the PCA clusters correctly by treatment along PC1.

Differential expression

We perform the differential expression analysis using DESeq2, shrinking the log2FoldChanges.

#Import the count table into a DESEq object
ddsmat_2_all = DESeqDataSetFromMatrix(counttable_2[,2:12], colData = coldata_2, design = ~ treatment)

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

#Estimate size factors: necessary to generate normalized counttable
ddsmat_2_all = estimateSizeFactors(ddsmat_2_all)

#Remove genes with 0 counttable in all samples to speed computation up
ddsmat_2_all_ne = ddsmat_2_all[rowSums(counts(ddsmat_2_all)) > 1, ]

#Differential expression
dds_2 = DESeq(ddsmat_2_all_ne)

## using pre-existing size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

#Extract results: non shrunk (ns) log2(FC) to compare with shrunk ones
res_2_ns_IgG_X203 = results(dds_2, contrast = c("treatment", "IgG", "X203"), alpha = 0.05)
res_2_ns_IgG_Saline = results(dds_2, contrast = c("treatment", "IgG", "Saline"), alpha = 0.05)
res_2_ns_X203_Saline = results(dds_2, contrast = c("treatment", "X203", "Saline"), alpha = 0.05)

#FC shrinkage
res_2_X203_IgG = lfcShrink(dds_2, contrast = c("treatment", "X203", "IgG"))

## using 'normal' for LFC shrinkage, the Normal prior from Love et al (2014).
## additional priors are available via the 'type' argument, see ?lfcShrink for details

res_2_IgG_Saline = lfcShrink(dds_2, contrast = c("treatment", "IgG", "Saline"))

## using 'normal' for LFC shrinkage, the Normal prior from Love et al (2014).
## additional priors are available via the 'type' argument, see ?lfcShrink for details

res_2_X203_Saline = lfcShrink(dds_2, contrast = c("treatment", "X203", "Saline"))

## using 'normal' for LFC shrinkage, the Normal prior from Love et al (2014).
## additional priors are available via the 'type' argument, see ?lfcShrink for details

#Create a results list
reslist_2 = list("res_2_X203_IgG" = res_2_X203_IgG, "res_2_IgG_Saline" = res_2_IgG_Saline, "res_2_X203_Saline" = res_2_X203_Saline)

for (i in names(reslist_2))
{
    reslist_2[[i]] = as.data.frame(reslist_2[[i]]) #conversion of DESeq results object to data frame
    reslist_2[[i]]$gene = mapping2[rownames(reslist_2[[i]]), 2] #add MGI gene symbols 
}

As a routine check, we look at the histogram of p values for each comparison, expecting to see the usual distribution (flat with a peak close to 0):

layout(matrix(1:4, nrow = 2))
for(i in names(reslist_2))
{
    hist(reslist_2[[i]]$pvalue, breaks = 100, col = "cadetblue", border = NA, main = paste("p value distribution for", i))
}

Distributions are as expected. We can now look at the differential expression profiles in volcano plots:

plot_volcano(reslist_2[[1]], alpha = 0.05, lfc = 1, maxdot = 40, cextest = 0.7, main = "X203 vs IgG", grid = F)

## Loading required package: plotrix

## 
## Attaching package: 'plotrix'

## The following object is masked from 'package:scales':
## 
##     rescale

## The following object is masked from 'package:gplots':
## 
##     plotCI

plot_volcano(reslist_2[[2]], alpha = 0.05, lfc = 1, maxdot = 40, cextest = 0.7, main = "IgG vs saline", grid = F)

plot_volcano(reslist_2[[3]], alpha = 0.05, lfc = 1, maxdot = 40, cextest = 0.7, main = "X203 vs saline", grid = F)

Gene Set Enrichment Analysis

Hallmark genesets

We begin by doing a GSEA using the Hallmark gene set from MSigDB. The input for the enrichment will be the stat column of the DESeq2 results object. We use 10^5 iterations and run the enrichment in parallel.

reslist_2_2 = list()
for (i in names(reslist_2))
{
    reslist_2_2[[i]] = reslist_2[[i]]
    reslist_2_2[[i]]$ensid = rownames(reslist_2_2[[i]])
    reslist_2_2[[i]] <- inner_join(reslist_2_2[[i]], bm, by=c("ensid"="ensembl_gene_id"))
}

for (i in names(reslist_2))
{
    write.table(file = paste("DEseq_2_", i, ".txt", sep = ""), reslist_2[[i]], sep = "\t", quote = F)
}
reslist_2 = lapply(reslist_2, function(x) x = x[-which(is.na(x$padj)),])


pathways = MSigDB$HALLMARK #select the Hallmark collection

fgsea_hallmark_list_2 = list() #create empty list to store GSEA results

for (i in 1:length(reslist_2))
{
res2 <- reslist_2_2[[i]] %>% 
  dplyr::select(hsapiens_homolog_associated_gene_name, stat) %>% 
  na.omit() %>% 
  distinct() %>% 
  group_by(hsapiens_homolog_associated_gene_name) %>% 
  summarize(stat=mean(stat))
ranks <- deframe(res2)

fgsea_hallmark_list_2[[names(reslist_2_2)[i]]] <- fgsea(pathways=pathways, stats=ranks, nperm = 100000, BPPARAM = SnowParam()) #do the enrichment analysis
fgsea_hallmark_list_2[[names(reslist_2_2[i])]] <- as.data.frame(fgsea_hallmark_list_2[[names(reslist_2_2[i])]][,c(1,3,5)])
rownames(fgsea_hallmark_list_2[[names(reslist_2_2)[i]]]) = fgsea_hallmark_list_2[[names(reslist_2_2[i])]]$pathway
fgsea_hallmark_list_2[[names(reslist_2_2[i])]]$pathway = NULL

}

## Warning in fgsea(pathways = pathways, stats = ranks, nperm = 1e+05, BPPARAM = SnowParam()): There are ties in the preranked stats (6.54% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

## Warning in fgsea(pathways = pathways, stats = ranks, nperm = 1e+05, BPPARAM = SnowParam()): There are ties in the preranked stats (6.49% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

## Warning in fgsea(pathways = pathways, stats = ranks, nperm = 1e+05, BPPARAM = SnowParam()): There are ties in the preranked stats (6.41% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

hallmark_2_df = as.data.frame(do.call(cbind, fgsea_hallmark_list_2))

hallmark_2_df_p = hallmark_2_df[,grep("padj", colnames(hallmark_2_df))]
hallmark_2_df_n = hallmark_2_df[,grep("NES", colnames(hallmark_2_df))]

par(mar = c(6, 1, 1, 20))
bubbleMap(valuedf = hallmark_2_df_n, pvaluedf = hallmark_2_df_p, color_pal = viridis(25, option = "A"), maplabel = "NES", cex.binned = F)

## Using func as id variables

Interestingly, we can detect the anti-inflammatory effect of the antibody by seeing how the inflammatory signatures are mostly down-regulated (significant, negative normalized enrichment score).

Gene Ontology genesets

We then use the whole GO collection to perform a similar analysis. However, since showing all the results would be impossible, we pre-select some interesting signatures (ECM and inflammation) and then we show, unbiasedly, the top 20 most up- and down-regulated signatures in each comparison.

pathways = MSigDB$C5_GENE_ONTOLOGY #select the Hallmark collection

fgsea_2_go_list = list() #create empty list to store GSEA results

for (i in 1:length(reslist_2_2))
{
res2 <- reslist_2_2[[i]] %>% 
  dplyr::select(hsapiens_homolog_associated_gene_name, stat) %>% 
  na.omit() %>% 
  distinct() %>% 
  group_by(hsapiens_homolog_associated_gene_name) %>% 
  summarize(stat=mean(stat))
ranks <- deframe(res2)

fgsea_2_go_list[[names(reslist_2_2)[i]]] <- fgsea(pathways=pathways, stats=ranks, nperm = 100000, BPPARAM = SnowParam()) #do the enrichment analysis
fgsea_2_go_list[[names(reslist_2_2[i])]] <- as.data.frame(fgsea_2_go_list[[names(reslist_2_2[i])]][,c(1,3,5)])
rownames(fgsea_2_go_list[[names(reslist_2_2)[i]]]) = fgsea_2_go_list[[names(reslist_2_2[i])]]$pathway
fgsea_2_go_list[[names(reslist_2_2[i])]]$pathway = NULL

}

## Warning in fgsea(pathways = pathways, stats = ranks, nperm = 1e+05, BPPARAM = SnowParam()): There are ties in the preranked stats (6.54% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

## Warning in fgsea(pathways = pathways, stats = ranks, nperm = 1e+05, BPPARAM = SnowParam()): There are ties in the preranked stats (6.49% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

## Warning in fgsea(pathways = pathways, stats = ranks, nperm = 1e+05, BPPARAM = SnowParam()): There are ties in the preranked stats (6.41% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

common_gsea_2 = intersect(intersect(rownames(fgsea_2_go_list[[1]]), rownames(fgsea_2_go_list[[2]])), rownames(fgsea_2_go_list[[3]]))
fgsea_2_go_list = lapply(fgsea_2_go_list, function(x) x[common_gsea_2,])
go_df_2 = as.data.frame(do.call(cbind, fgsea_2_go_list))

go_df_p_2 = go_df_2[,grep("padj", colnames(go_df_2))]
go_df_n_2 = go_df_2[,grep("NES", colnames(go_df_2))]

inflammation_rows_2 = grep("INFLAM", rownames(go_df_2))
ecm_rows = grep("EXTRACELLULAR_MATRIX", rownames(go_df_2))

par(mar = c(6, 1, 1, 20))

bubbleMap(valuedf = go_df_n_2[c(ecm_rows,inflammation_rows_2),], pvaluedf = go_df_p_2[c(ecm_rows,inflammation_rows_2),], maplabel = "NES", color_pal = viridis(25, option = "A"), cex.binned = F)

## Using func as id variables

Top 20 most up- and down-regulated lollipop plot in each comparison:

tdf = as.data.frame(cbind(go_df_n_2[,1], go_df_p_2[,1]))
rownames(tdf) = rownames(go_df_n_2)
colnames(tdf) = c("NES", "padj")
tdf$pathway = rownames(tdf)
lollipop_GSEA(tdf, show.max = 20, main = "X203 vs IgG GSEA", xlim = c(-5, 5))

tdf = as.data.frame(cbind(go_df_n_2[,2], go_df_p_2[,2]))
rownames(tdf) = rownames(go_df_n_2)
colnames(tdf) = c("NES", "padj")
tdf$pathway = rownames(tdf)
lollipop_GSEA(tdf, show.max = 20, main = "IgG vs Saline GSEA", xlim = c(-5, 5))

tdf = as.data.frame(cbind(go_df_n_2[,3], go_df_p_2[,3]))
rownames(tdf) = rownames(go_df_n_2)
colnames(tdf) = c("NES", "padj")
tdf$pathway = rownames(tdf)
lollipop_GSEA(tdf, show.max = 20, main = "X203 vs Saline GSEA", xlim = c(-5, 5))

Specific ECM regulation by X203

We can show that both IgG and X203 target more ECM than inflammation by plotting the normalized enrichment score of X203 vs IgG against the score of IgG vs Saline. The fact that magenta dots are more to the upper left corner of the plot thatn the orange ones means that ECM is more affected than inflammation.

plot( go_df_n_2[c(ecm_rows,inflammation_rows_2),1], go_df_n_2[c(ecm_rows,inflammation_rows_2),2], pch = 16, col = c(rep("darkmagenta", length(ecm_rows)), rep("orange", length(inflammation_rows_2))), ylim = c(-3,3), xlim = c(-3,3), xlab = "X203 vs IgG NES", ylab = "IgG vs Saline NES", main = "Enrichment score comparison")
abline(h = 0) + abline(v = 0) + abline(0,1, lty = 2) + abline(0,-1, lty = 2)

## integer(0)

legend ("topright", bty = "n", pch = 16, col = c("darkmagenta", "orange"), legend = c("ECM", "Inflammation"))

As another confirmation, we can plot the distribution of log2FoldChanges in X203 and IgG of Gene Ontology ECM and inflammation genes (dots show significant genes):

go_inf = MSigDB$C5_GENE_ONTOLOGY$GO_INFLAMMATORY_RESPONSE
go_inf_mouse = bm[bm[,2] %in% go_inf,1]

lfc_X203_IgG_inf = reslist_2$res_2_X203_IgG[go_inf_mouse,2]
lfc_X203_IgG_inf_sig = reslist_2$res_2_X203_IgG[intersect(go_inf_mouse, rownames(reslist_2$res_2_X203_IgG[which(reslist_2$res_2_X203_IgG$padj < 0.05),])),2]

go_ecm = MSigDB$C5_GENE_ONTOLOGY$GO_EXTRACELLULAR_MATRIX_COMPONENT
go_ecm_mouse = bm[bm[,2] %in% go_ecm,1]

lfc_X203_IgG_ecm = reslist_2$res_2_X203_IgG[go_ecm_mouse,2]
lfc_X203_IgG_ecm_sig =  reslist_2$res_2_X203_IgG[intersect(go_ecm_mouse, rownames(reslist_2$res_2_X203_IgG[which(reslist_2$res_2_X203_IgG$padj < 0.05),])),2]


boxplot(lfc_X203_IgG_ecm, lfc_X203_IgG_inf, ylab = "log2FoldChange", border = c("purple", "orange"), xaxt = "n", ylim =c(-2,2))
axis(1, at = 1:2, labels = c("ECM", "INF"))
beeswarm(lfc_X203_IgG_ecm_sig, at = 1, add = T, pch = 21, bg = "purple")
beeswarm(lfc_X203_IgG_inf_sig, at = 2, add = T, pch = 21, bg = "orange")
abline(h = 0, lty = 2)

And their volcano plots:

layout(matrix(1:2, nrow = 1))
plot_volcano(reslist_2$res_2_X203_IgG[go_ecm_mouse,], grid =F, cextest = 0.7, main = "ECM genes (X203 vs IgG)", xlim = c(-2,2))
plot_volcano(reslist_2$res_2_X203_IgG[go_inf_mouse,], grid =F, cextest = 0.7, main = "Inflammation genes (X203 vs IgG)", xlim = c(-2,2))

Taken together, these results show that the effect on ECM is more prominent than that on inflammation, although the number of significant differentially expressed genes (at an FDR threshold of 0.05) is similar. The enrichment depends on the size of the respective categories. However, the distribution of all fold changes is more shifted to a negative median in ECM genes compared to inflammation genes.

To better represent this result we isolate from the GSEA in X203 vs IgG the inflammation and ECM genesets and we show the lollipop plot:

tdf = as.data.frame(cbind(go_df_n_2[,1], go_df_p_2[,1]))
rownames(tdf) = rownames(go_df_n_2)
colnames(tdf) = c("NES", "padj")
tdf$pathway = rownames(tdf)
tdf = tdf[c(inflammation_rows_2, ecm_rows),]
par(mar = c(5,1,3,1))
lollipop_GSEA(tdf, show.max = 28, main = "ECM and Inflammation X203 vs IgG GSEA", xlim = c(-5, 5))

We can see how the only significant inflammation categories have a lower enrichment than the significant ECM categories.

Session info

sessionInfo()

## R version 3.5.2 (2018-12-20)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.4
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] C/UTF-8/C/C/C/C
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] plotrix_3.7-5               fgsea_1.8.0                
##  [3] Rcpp_1.0.1                  MSigDB_1.0.6.1             
##  [5] colorspace_1.4-1            scales_1.0.0               
##  [7] treemapify_2.5.3            viridis_0.5.1              
##  [9] viridisLite_0.3.0           pheatmap_1.0.12            
## [11] gplots_3.0.1.1              beeswarm_0.2.3             
## [13] biomaRt_2.38.0              forcats_0.4.0              
## [15] stringr_1.4.0               dplyr_0.8.0.1              
## [17] purrr_0.3.2                 readr_1.3.1                
## [19] tidyr_0.8.3                 tibble_2.1.1               
## [21] ggplot2_3.1.1               tidyverse_1.2.1            
## [23] edgeR_3.24.3                limma_3.38.3               
## [25] sva_3.30.1                  genefilter_1.64.0          
## [27] mgcv_1.8-28                 nlme_3.1-139               
## [29] DESeq2_1.22.2               SummarizedExperiment_1.12.0
## [31] DelayedArray_0.8.0          BiocParallel_1.16.6        
## [33] matrixStats_0.54.0          Biobase_2.42.0             
## [35] GenomicRanges_1.34.0        GenomeInfoDb_1.18.2        
## [37] IRanges_2.16.0              S4Vectors_0.20.1           
## [39] BiocGenerics_0.28.0        
## 
## loaded via a namespace (and not attached):
##  [1] htmlTable_1.13.1       XVector_0.22.0         base64enc_0.1-3       
##  [4] rstudioapi_0.10        ggfittext_0.6.0        bit64_0.9-7           
##  [7] AnnotationDbi_1.44.0   lubridate_1.7.4        xml2_1.2.0            
## [10] splines_3.5.2          geneplotter_1.60.0     knitr_1.22            
## [13] Formula_1.2-3          jsonlite_1.6           broom_0.5.2           
## [16] annotate_1.60.1        cluster_2.0.8          compiler_3.5.2        
## [19] httr_1.4.0             backports_1.1.4        assertthat_0.2.1      
## [22] Matrix_1.2-17          lazyeval_0.2.2         cli_1.1.0             
## [25] acepack_1.4.1          htmltools_0.3.6        prettyunits_1.0.2     
## [28] tools_3.5.2            gtable_0.3.0           glue_1.3.1            
## [31] GenomeInfoDbData_1.2.0 fastmatch_1.1-0        cellranger_1.1.0      
## [34] gdata_2.18.0           xfun_0.6               rvest_0.3.3           
## [37] gtools_3.8.1           XML_3.98-1.19          zlibbioc_1.28.0       
## [40] hms_0.4.2              RColorBrewer_1.1-2     curl_3.3              
## [43] yaml_2.2.0             memoise_1.1.0          gridExtra_2.3         
## [46] rpart_4.1-15           reshape_0.8.8          latticeExtra_0.6-28   
## [49] stringi_1.4.3          RSQLite_2.1.1          checkmate_1.9.1       
## [52] caTools_1.17.1.2       rlang_0.3.4            pkgconfig_2.0.2       
## [55] bitops_1.0-6           evaluate_0.13          lattice_0.20-38       
## [58] htmlwidgets_1.3        bit_1.1-14             tidyselect_0.2.5      
## [61] plyr_1.8.4             magrittr_1.5           R6_2.4.0              
## [64] snow_0.4-3             generics_0.0.2         Hmisc_4.2-0           
## [67] DBI_1.0.0              pillar_1.3.1           haven_2.1.0           
## [70] foreign_0.8-71         withr_2.1.2            survival_2.44-1.1     
## [73] RCurl_1.95-4.12        nnet_7.3-12            modelr_0.1.4          
## [76] crayon_1.3.4           KernSmooth_2.23-15     rmarkdown_1.12        
## [79] progress_1.2.0         locfit_1.5-9.1         grid_3.5.2            
## [82] readxl_1.3.1           data.table_1.12.2      blob_1.1.1            
## [85] digest_0.6.18          xtable_1.8-3           munsell_0.5.0